Monte-Carlo simulations of the clean and disordered contact process in three 

dimensions 
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The absorbing-state transition in the three-dimensional contact process with and without 
quenched randomness is investigated by means of Monte-Carlo simulations. In the clean case, 
a reweighting technique is combined with a careful extrapolation of the data to infinite time to 
determine with high accuracy the critical behavior in the three-dimensional directed percolation 
universality class. In the presence of quenched spatial disorder, our data demonstrate that the 
absorbing-state transition is governed by an unconventional infinite-randomness critical point fea- 
turing activated dynamical scaling. The critical behavior of this transition does not depend on 
the disorder strength, i.e., it is universal. Close to the disordered critical point, the dynamics is 
characterized by the nonuniversal power laws typical of a Griffiths phase. We compare our findings 
to the results of other numerical methods, and we relate them to a general classification of phase 
transitions in disordered systems based on the rare region dimensionality. 

PACS numbers: 05.70.Ln, 64.60.Ht, 02.50.Ey 



I. INTRODUCTION 

The macroscopic behavior of many-particle systems 
far from equilibrium can abruptly change when an ex- 
ternal parameter is changed. The resulting nonequilib- 
rium phase transitions separate different nonequilibrium 
steady states. They are characterized by strong fluctu- 
ations and cooperative phenomena over large distances 
and times, analogous to the behavior at equilibrium 
phase transitions. Examples of nonequilibrium phase 
transitions can be found in catalytic reactions, growing 
interfaces, turbulence, and traffic jams as well as in the 
dynamics of epidemics and other biological populations 
(see, e.g., Refs. [JjjJ). 

A well-studied class of nonequilibrium phase transi- 
tions are the absorbing-state transitions between active, 
fluctuating steady states and inactive, absorbing states 
in which fluctuations cease completely. The generic uni- 
versality class for absorbing-state transitions is the di- 
rected percolation (DP) class Janssen and Grass- 
berger [ToL fllj conjectured that all absorbing-state tran- 
sitions with a scalar order parameter and short-range in- 
teractions belong to this class, provided they do not fea- 
ture extra symmetries or conservation laws. Additional 
symmetries or conservation laws can lead to other uni- 
versality classes such as the parity conserving class or 
^-symmetric directed percolation (DP2) (see, e.g., Refs. 

a ft- 

Although absorbing state transitions are ubiquitous in 
theory and computer simulations, experimental observa- 
tions of their universality classes were lacking for a long 
time [12j. A full verification of the DP universality class 
was recently achieved in the transition between two tur- 
bulent states in a liquid crystal (l3^ . Other absorbing 
state transitions were found in periodically driven sus- 
pensions [Til . [TBj and in superconducting vortices [lfjj . 

In many experimental systems, one can expect impuri- 
ties and defects to play an important role. Indeed, it has 



been suggested [H| that such quenched spatial disorder is 
one of the key reasons for the surprising rarity of the DP 
universality class in experiments. The influence of dis- 
order on absorbing state transitions is therefore a prime 
problem in the field. According to the Harris criterion 
|17| , a clean critical point is stable against the introduc- 
tion of weak spatial disorder if its correlation length crit- 
ical exponent v±_ fulfills the inequality dv±_ > 2 where d 
is the space dimensionality. The values of v±_ in the clean 
DP universality class are approximately 1.1 in one dimen- 
sion, 0.73 in two dimensions, and 0.58 in three dimensions 
[3| ■ The Harris criterion is thus violated, and spatial dis- 
order is expected to change the critical behavior. This 
heuristic result was confirmed by a field-theoretic renor- 
malization group study [l8[ which found runaway flow 
towards large disorder. Early Monte-Carlo simulations 
[1914251 ] demonstrated unusually slow dynamics but could 
not resolve the ultimate fate of the transition. 

In recent years, a comprehensive understanding of 
the one-dimensional disordered contact process has been 
achieved by a combination of analytical and numerical 
approaches. A stron g-d isorder renormalization group 
(SDRG) analysis [H, Wf\ established that the critical 
point is of exotic infinite-randomness type and charac- 
terized by activated (exponential) dynamical scaling. It 
belongs to the same universality class as the random 
transverse-field Ising chain 0, [29| , at least for suffi- 
ciently strong disorder. These predictions were confirmed 
by large-scale Monte-Carlo simulations [3(| that also pro- 
vided evidence for the critical behavior being universal, 
i.e., independent of the disorder strength. In higher di- 
mensions, the SDRG cannot be solved analytically. How- 
ever, by using a numerical implementation of the SDRG 
[3II ]. the infinite-randomness scenario was found to be 
valid in two dimensions, in agreement with Monte- Carlo 
simulations of the contact process on diluted lattices 

Here, we extend our Monte-Carlo simulations of the 
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contact process on diluted lattices to three space dimen- 
sions. Using large lattices of up to 999 3 sites and long 
times up to 10 s , we provide strong evidence for the non- 
equilibrium phase transition of the disordered contact 
process being governed by an infinite-randomness critical 
point. We determine the critical exponents and find them 
to be universal, i.e., independent of disorder strength. In 
contrast, the dynamics in the Griffiths region between 
the clean and disordered critical points is characterized 
by nonuniversal power laws. As a byproduct of our sim- 
ulations, we also obtain high-precision estimates for the 
critical exponents of the clean contact process in three 
dimensions. 

The paper is organized as follows. The contact pro- 
cess on a diluted lattice is introduced in Sec. HU We 
briefly summarize the scaling theories for conventional 
and infinite-randomness critical points in Sec. Mil In Sec. 
IIV1 we describe our simulation method and present the 
results. We conclude in Sec. [V] 



II. DEFINITION OF THE CONTACT PROCESS 

The contact process [35| can be viewed as a model for 
the spreading of an epidemic in space. Consider a hyper- 
cubic d-dimensional lattice of L d sites. Each site can be 
in one of two states, either active (infected) or inactive 
(healthy). The time evolution of the contact process is 
a continuous-time Markov process during which infected 
sites heal spontaneously at a rate fi while healthy sites 
become infected by their neighbors at a rate An/ (2d). 
Here, n is the number of sick nearest neighbors of the 
given site. The infection rate A and the healing rate \x 
(which can be set to unity) are the external control pa- 
rameters that govern the behavior of the system. 

The steady states of the contact process can be easily 
understood at a qualitative level. For A <C fx, healing 
dominates over infection, and the epidemic eventually 
dies out completely. The system therefore always ends up 
in the absorbing steady state without any infected sites. 
This is the inactive phase. In contrast, the density of 
infected sites remains nonzero in the long-time limit if the 
infection rate A is sufficiently large, i.e., the system is in 
the active phase. The nonequilibrium transition between 
these two phases, which occurs at a critical infection rate 
A°, belongs to the DP universality class. 

Quenched spatial disorder can be introduced into the 
contact process in different ways, e.g., by making the in- 
fection and healing rates random variables, or by using 
a random lattice instead of a regular one. Here, we ran- 
domly dilute the regular lattice by removing each site 
with probability p 36]. In the context of an epidemic, 
a vacancy can be interpreted as a site that is immune 
against the infection. For vacancy concentrations p be- 
low the percolation threshold p c , the lattice still has an 
infinite connected cluster of sites that can support an 
active phase of the contact process. If the vacancy con- 
centration is above p c , an infinite cluster does not ex- 
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FIG. 1. (Color online) Phase diagram of the contact process 
on a site-diluted cubic lattice (inverse critical infection rate 
A^ 1 vs vacancy concentration p). MCP marks the multicriti- 
cal point. The black dots show the actual simulation results, 
the lines are guides to the eye. 

ist. Instead, the lattice consists of disconnected finite-size 
clusters. As the epidemic dies out on any finite cluster 
in the long-time limit, an active phase is impossible for 
p > p c . This leads to the phase diagram shown in Fig. [1] 
Specifically, there are two different nonequilibrium phase 
transitions: (i) the so-called generic transition for p < p c , 
driven by the dynamic fluctuations of the contact pro- 
cess and (ii) the lattice percolation transition occurring 
at p = p c for sufficiently large infection rates [37l l38j] ■ The 
two phase transition lines meet at a multicritical point 

The central quantity of the contact process is the den- 
sity of infected sites at time t, 

K*) = 7^£K(*)> ■ (1) 

r 

Here, n r (t) is the occupation of site r at time t, i.e., 
n r (t) = I if the site is infected and n r (t) — if it is 
healthy. (...) denotes the average over all realizations 
of the Markov process. The order parameter of the 
absorbing-state phase transition is given by the steady 
state density 

Pstat = lim p(t) . (2) 



III. SCALING THEORIES OF ABSORBING 
STATE TRANSITIONS 

In this section, we summarize the scaling theories of the 
nonequilibrium transitions in the clean and disordered 
contact process to the extent necessary for analyzing our 
Monte-Carlo data. We contrast the cases of conventional 
power-law scaling and activated scaling. More details can 
be found, for instance, in Ref. [3] for the power-law case 
and in Refs. |4l| for the activated case. 
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A. Conventional critical points 

The DP universality class is characterized by three in- 
dependent critical exponents, for example, ft, and 
z. The order parameter exponent f3 describes how the 
steady state density varies as the infection rate A ap- 
proaches its critical value A c from above, 

Pstat ~ (A - \ c f ~ A' 3 . (3) 

Here, A = (A — A c )/A c is the dimensionless distance from 
criticality. The correlation length exponent v±_ describes 
the divergence of the correlation length £j_ at criticality, 

U ~ |A|-^ . (4) 

The correlation time £|| diverges like a power of the cor- 
relation length, 

(5) 

which defines the dynamical exponent z. In terms of 
these exponents, the scaling form of the density as a func- 
tion of A, time i, and system size L reads 

p(A, t, L) = b^^piAb- 1 /^ , tb z , Lb) . (6) 

Here, b is an arbitrary dimensionless length scale factor. 

If the time evolution starts at time from a single in- 
fected site in an otherwise inactive lattice, one can ask 
what is the probability that an active cluster survives at 
time t. In the DP universality class, this survival proba- 
bility P s has the same scaling form as the density [421 ] . 

P s (A,t,L) = bP! v± P s (M- l / v ^,tb z ,Lb) . (7) 

The correlation (or pair connectedness) function 
C(r, t) = {n T (t)n (0)} is given by the probability that 
site r is infected at time t when the time evolution starts 
from a single infected site at r = and time 0. The scale 
dimension of C is 2/3 /vj_ because it involves a product of 
two densities, leading to the scaling form [431 ] 

C{A 1 r,t,L) = b 2 ^^C{Ab~ 1 /^,rb,tb z ,Lb) . (8) 

The total number N s of sites in the active cluster can be 
calculated by integrating the correlation function over all 
space, resulting in 

N S (A, t, L) = b 2p/ ^- d N s {Ab- 1/l '^ , tb z ,Lb) . (9) 

The mean-square radius R of the active cluster has the 
dimension of a length. Its scaling form therefore reads 

R(A, t, L) = 6- 1 i?(Afe- 1 ^ , tb\ Lb) . (10) 

The functional dependencies of p, P s , N s and R on the 
parameters A, i, and L can be easily derived from the 
scaling forms by setting the scale factor b to appropriate 
values. This leads to the following time dependencies at 
the critical point A = and in the thermodynamic limit 



L — >• oo. In the long-time limit, the density of infected 
sites and the survival probability obey the power laws 

p(t)~t~ s , P s (t)~t~ s (11) 

with S = /3/(v±z). The mean-square radius and number 
of infected sites of a cluster starting from a single seed 
site behave as 

Rtf-t 1 /*, N s (t)~t e (12) 

where = d/z — 2/3 /(v±z) is the so-called critical initial 
slip exponent. By taking the derivative of eqs. ©, (|7|), 
and |9]) with respect to A, we also find that 

dlnp d\nP s d\nN s 1/(u±z) , , 
OA ~ dA ~ OA 1 J 

which will be useful for measuring v±. 

B. Infinite-randomness critical points 

Infinite-randomness critical points feature extremely 
slow dynamics, represented by an exponential (activated) 
relation between correlation length and time 

ln(£||Ao)~& (14) 

rather than the power-law dependence ([S]). It is charac- 
terized by the so-called tunneling exponent ip, and to is 
a nonuniversal microscopic time scale. The exponential 
relation between time and length scales implies that the 
dynamical exponent z is formally infinite. In contrast to 
the dynamical scaling, the static scaling relations remain 
of power-law type, i.e., eqs. ([3]) and (0]) remain valid. 

The scaling forms of disorder-averaged observables can 
be obtained by simply substituting the variable combi- 
nation ln(t/to)& for tb z in the arguments of the scaling 
functions: 

p(AM(t/t Q ),L)^b^^p(Ab- 1 /^Mt/toP',Lb) , 

(15) 

P s {AMt/to),L) = b /^P s (Ab- 1 ^Mt/to)b^,Lb) , 

(16) 

N(A,]n(t/to),L) = b 2 ^ /v± - d N(Ab- 1/v± ,]n{t/to)b^,Lb) , 

(17) 

R{A,bx(t/to),L) =b- 1 R(Ab- 1/l,± .\n(t/t )b' p ,Lb) . 

(18) 

Consequently, the critical time dependencies of the 
density of active sites and the survival probability (in 
the thermodynamic limit) are logarithmic, 

p(t) ~ [ln(t/t )]- S , P s (t) ~ [Ht/t )]- S (19) 

with (5 = (3/(vj_iJj). The radius and number of active sites 
in a cluster starting from a single seed site vary as 

R(t) ~ [Ht/to)} 1 ^, N s (t) ~ [\n(t/t )f (20) 
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with 6 = d/ip — 2f3/(i'±ip). Taking the derivatives of eqs. 
([15]). (|T5j). and dnjl with respect to A yields 



dlnp d\nP s d\nN s 



dA 



dA 



dA 



Mt/to)] 1 /^ . (21) 



C. Griffiths region 

In the presence of spatial disorder, the contact process 
displays unconventional behavior not just at the criti- 
cal point but also in its vicinity because rare active spa- 
tial regions dominate the long-time dynamics. This phe- 
nomenon is an example of the well-known Griffiths sin- 
gularities [IH that generally occur at phase transitions 
in disordered systems (see Ref. [4(3] for a review). The 
Griffiths singularities in the spatially disordered contact 
process can be understood as follows [l9j . 

The inactive phase must be divided into two regions, 
(i) If the infection rate is below the clean critical value, 
A < A°, the behavior is conventional. This means, the 
system approaches the absorbing state exponentially fast. 
The decay time increases with A and diverges as |A — 
;\0|— »v± where z and v±_ are the exponents of the clean 
critical point [3(1 5H . 

(ii) If the infection rate is in the so-called Griffiths 
region (or Griffiths phase) between the clean and dirty 
critical values, A° < A < A c , the system is globally still 
in the inactive phase (i.e., it eventually decays into the 
absorbing state). However, in the thermodynamic limit, 
one can find arbitrarily large spatial regions devoid of 
vacancies. These rare regions are locally in the active 
phase. Although they cannot support a non-zero steady 
state density because they are of finite size, their time 
decay is very slow as it requires a rare, exceptionally 
large density fluctuation. 

The contribution of the rare regions to the density of 
infected sites can be expressed as the integral 



Pit) 



dL r Lf w(L r )exp[—t/r(L r )] 



(22) 



Here, w denotes the probability for finding a spatial re- 
gion of size L r that does not contain any vacancies, and 
r(L r ) is the life time of the contact process on such a 
rare region. Basic combinatorics gives 



w{L r ) ~ exp(-pL^) 



(23) 



with p = — ln(l — p) (up to pre-exponential factors). In 
the Griffiths phase, the life time of a rare region depends 
exponentially on its volume, 



r(L r ) - exp(aLf) 



(24) 



because a coordinated fluctuation of the entire region is 
necessary to take it to the absorbing state [ill EH IH[ . 
The constant a vanishes at the clean critical infection rate 
A° and increases with A. Evaluating the integral (|2"2"|) in 
saddle-point approximation, we obtain a power-law time 



dependence for the density. The survival probability P s 
shows exactly the same time dependence, 



—p/a j — djz 



(25) 



where z' = da/p is the nonuniversal dynamical expo- 
nent in the Griffiths region. The behavior of z' close to 
the dirty critical point can be obtained from the SDRG 
analysis [13, [H, [3l| . As A approaches A c , z' diverges as 



z'~|A-A c |-^ 



(26) 



where tp and u± are the critical exponents of the infinite- 
randomness critical point. 



IV. MONTE-CARLO SIMULATIONS 
A. Simulation method 

To perform Monte Carlo simulations of the contact 
process on randomly diluted cubic lattices, we followed 
the implementation described, for instance, by Dickman 
47| . The algorithm starts at time t — from some con- 
figuration of infected and healthy sites and consists of a 
sequence of events. During each event an infected site 
is randomly chosen from a list of all N a infected sites, 
then a process is selected, either infection of a neighbor 
with probability A/(l + A) or healing with probability 
1/(1 + A). For infection, one of the six neighbor sites is 
chosen at random. The infection succeeds if this neigh- 
bor is healthy (and not a vacancy site). The time is then 
incremented by 1/N a . 

Using this algorithm, we simulated systems with sizes 
of up to 999 3 sites and vacancy concentrations p = 0, 0.2, 
0.3, 0.4, 0.5, 0.6 and p c = 0.6883920 48]. To cope with 
the slow dynamics of the disordered contact process, we 
simulated long times up to 10 8 . All results were averaged 
over a large number of disorder configurations, precise 
numbers will be given below. 

We carried out two different types of simulations, (i) 
The majority of runs started from a single active site in 
an otherwise inactive lattice (spreading runs); we moni- 
tored the survival probability P s (t) , the number of sites 
N s (t) of the active cluster, and its radius R(t). (ii) For 
comparison, we also performed a few runs that started 
from a completely active lattice during which we observed 
the time evolution of the density p(t). 

We employed two different high-quality, long-period 
random number generators. Most simulations used 
LFSR113 proposed by L'Ecuyer [49j]. We verified the 
validity of the results by means of the 2005 version of 
Marsaglia's KISS [5(| ■ The total computational effort 
for the work described in this paper was about 100 000 
CPU days on the Pegasus cluster at Missouri S&T. 

Figure [T] gives an overview of the phase diagram re- 
sulting from these simulations. As expected, the critical 
infection rate A c increases with increasing impurity con- 
centration. 
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FIG. 2. (Color online) Effective critical exponent —<d/S — 
(dlnN s )/(dlnP a ) vs. <" 1/2 . The critical curve is marked by 
dots, and the dashed line is a linear extrapolation to t = oo. 



B. Contact process on an undiluted lattice 

The purpose of studying the clean three-dimensional 
contact process is two-fold, (i) to test our implementa- 
tion of the contact process and (ii) to compute highly 
accurate estimates of the critical exponents in the three- 
dimensional DP universality class. To reduce the numer- 
ical effort, we applied the clever reweighting technique 
proposed in Ref. [47j |. 

After a few test calculations aimed at bracketing the 
critical point, we performed two large spreading runs 
(starting from a single active site) at A = 1.3168400. By 
reweighting with a step AA = 0.0000025, we generated 
data for A between 1.3168150 and 1.3168650. The first 
run consisted of 4 x 10 8 trials using the LFSR113 random 
number generator, the second consisted of 5 x 10 8 trials 
using the KISS random number generator. The maxi- 
mum time of both runs was 5 x 10 4 . As the data of both 
runs agree within their statistical error, we averaged their 
results. The system size, 850 3 sites, was chosen such that 
the active cluster stayed smaller than the sample during 
the entire time evolution, eliminating finite-size effects. 

To find the location of the critical point and to measure 
the critical exponents, we define effective (running) ex- 
ponents via the logarithmic derivatives of various observ- 
ables. These effective exponents are then extrapolated to 
t = oo. The finite-size scaling exponent /3/V-L (scale di- 
mension of the order parameter) can be determined from 
the relation between N s and P s . Combining (fTTj) and 
(TT2]) yields N s ~ P~ e/S with 0/5 = 3v ± /fi - 2. Fig- 
ure [5] shows the effective exponent (dlnAT s )/(dlnP s ) as 
a function of t~ v with y = 1/2. (The value 1/2 was cho- 
sen empirically to allow a linear extrapolation to t = oo.) 
From this plot, we estimate the critical infection rate to 
be 
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FIG. 3. (Color online) Effective critical exponents 
l/z = (din R)/ (dint), 5 = (din P s ) /(dint) and 8 = 
(din N a )/ (dint) vs. t~ 1/2 . The critical curves are marked by 
dots, and the dashed lines are linear extrapolations to t — oo. 



We verified this value by performing an extra run directly 
at A = 1.316835 using 4 x 10 8 trials on a system of size 
999 3 with a maximum time of 10 5 . 

Extrapolating the effective exponent to t = oo, we ob- 
tain 0/8 = 0.1442(3) sys (2) ran where the values in brack- 
ets represent estimates of the systematic and random er- 
rors of the last digit. The systematic error stems from 
the uncertainties of AJ? and the extrapolation exponent y 
while the random error is due to the Monte-Carlo noise. 
The resulting value of the finite-size scaling exponent is 
/3/fX = 1-3991(4). An estimate for this exponent can 
also be obtained from the relation between iV s and R. 
Extrapolating the effective exponent as above yields the 
identical value f}/v± = 1.3991(4). 

To determine the exponents z, 8, and 0, we apply the 
same type of analysis to the logarithmic derivatives of 
R, P s , and N s with respect to time. The corresponding 
graphs are shown in Fig.[3J Extrapolation to t = oo yields 
the dynamical exponent l/z = 0.5267(l) sys (l) ran as well 
as S = 0.7367(5) sys (l) ran and 9 = 0.1062(2) sys (2) ran . 
These values fulfill hyperscaling because + 26 — 3/z — 
—0.0005(22), in excellent agreement with the exact result 
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Value This work Ref. [51] Ref. [47] Ref. [52] Ref. [53] 





1.316835(1) 


(*) 


1.31686(1) 1.31683(2) 1.3168(1) 


/3/V x 1.3991(4) 


1.395(4) 


1.394(1) 


1.392(5) 


v± 


0.5826(9) 




0.580(3) 


0.584(6) 


P 


0.815(2) 




0.808(5) 0.78(1) 


0.813(11) 


6 


0.7367(6) 




0.7263(11) 


0.732(4) 


e 


0.1062(4) 




0.110(1) 


0.114(4) 


z 


1.8986(8) 


1.916(5) 


1.919(4) 


1.901(5) 


2/2 


1.0534(4) 




1.042(2) 


1.052(3) 


viz 


1.106(2) 




1.114(4) 


1.11(1) 


0,2 


0.2016(6) 


0.216(3) 






£>/ 


1.6009(4) 




1.56(3) 





TABLE I. Critical infection rate and critical exponents of the 
clean three-dimensional contact process. The upright num- 
bers are directly given in the respective papers, the italic ones 
were calculated using scaling relations. The fractal dimension 
D f = 3 - fi/v x . (*)The authors of Ref. [HI used the value 
of A c found in Ref. [47[ as an input. 



of zero. 

Finally, we measure the exponent combination l/(vj_z) 
by analyzing the time dependencies of (d In P S )J (d\) and 
(d\n N s )/(dX) according to fll"3l) . Extrapolating the ef- 
fective exponent to t — oo as above yields l/{v±z) = 
0.9040(5) sys (5) ran . The correlation length and order pa- 
rameter exponents can be calculated by combining this 
value with our results for z and P/v± yielding v± = 
0.5826(9) and /3 = 0.815(2). 

In TablelH we compare our estimates for the critical ex- 
ponents with earlier results. The present estimates have 
significantly higher precision than the values in the lit- 
erature. They are roughly compatible with Jensen's val- 
ues [53[ within their given errors (for O, the difference 
is about twice the given error, though). However, they 
are clearly not compatible with the values given in Refs. 
[47l ] and [5l| (for 5, the difference is about ten times the 
given error, and for z it is about five times the given er- 
ror). We believe, this discrepancy can be traced back to 
the location of the critical point. According to our data, 
the infection rate A = 1.31686(1), identified as critical in 
Ref. [13] and also employed in [El]], is on the active side 
of the transition (it differs from our estimate by about 
three times the given error). As the survival probability 
decays more slowly in the active phase than at critical- 
ity, this may be responsible for the low 5- value and, via 
the hyperscaling relation, for the high z-value reported 
in Ref. S3. 



C. Contact process on a diluted lattice 

The remainder of Sec. IIVI focuses on the contact pro- 
cess on a diluted lattice. We tried to use the same 
reweighting technique as in the clean case to save com- 




FIG. 4. (Color online) Survival probability P s and number of 
active sites N s vs. time t for impurity concentration p = 0.5 
and several infection rates A. The critical curve at A c = 2.6906 
is marked by dots. 



puter time. However, these attempts were not success- 
ful. The reweighing method of Ref. [I?) considers a set 
of simulation runs (particular realizations of the Markov 
process) at some infection rate A and reweighs their sta- 
tistical probabilities according to a slightly different A'. 
This only works as long as the two infection rates are suf- 
ficiently close such that their sets of possible runs overlap 
significantly. This overlap decreases with increasing sim- 
ulation time. In the presence of disorder, particularly 
long simulation times are required because the critical 
dynamics is logarithmically slow. Thus, reweighting is 
restricted to very narrow A-intervals (too narrow com- 
pared to the range of infection rates we needed to ex- 
plore to determine the critical point). All results were 
thus obtained in the conventional manner by performing 
a separate run for each A- value. 

Figure [4] gives an overview over spreading simulations 
(starting from a single active seed site) for a vacancy con- 
centration p — 0.5. The data represent averages over at 
least 5000 disorder configurations, with 128 trials start- 
ing from random seed sites for each configuration. A 
system size of 500 3 sites ensured that the active cluster 
stayed smaller than the sample for the entire simulation 
run. The figure shows that the dynamics in the vicinity 
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FIG. 5. (Color online) N s versus P a for vacancy concentration 
p = 0.5 and several infection rates A close to the critical point. 



of the phase transition is very slow. In particular, the 
time-dependence of the survival probability appears to 
be slower than a power law, in agreement with the acti- 
vated scaling scenario of Sec. IIIIBI Moreover, the data 
show indications of Griffiths singularities, i.e., nonuni- 
versal power-law behavior somewhat below the critical 
infection rate. We also note that the number of sites in 
the active cluster N s decreases with time at the transi- 
tion, in contrast to the clean case and to the diluted case 
in two dimensions [33[ • This implies a negative exponent 

e. 

To find the precise location of the critical point within 
the activated scaling scenario, one might be tempted to 
search for power-law relations between lni and observ- 
ables such P s and N s (either by plotting In P s and In N s 
vs. In hit or by analyzing the corresponding effective ex- 
ponents). However, this method is highly unreliable as 
the unknown microscopic time scale to in (fH?)) and ([20]) 
provides a correction to scaling via ln(t/to) = Int — Into- 
This strongly influences the results because the simula- 
tions cover only a moderate r ang e in Inf. (In the two- 
dimensional simulations, Ref. [33[, it was found that ne- 
glecting to could change the apparent value of S from its 
correct value of 1.9 to 3.) 

To circumvent this problem, we follow the method de- 
vised in Ref. [33|]. It is based on the observation that 
to has the same value in the scaling forms of all quan- 
tities because it is related to the energy scale Qq of the 
underlying SDRG. Consequently, if one analyzes the re- 
lation between N s and P s or other such combinations of 
observables, the critical point corresponds to power-law 
behavior (independent of the value of to) as long as all 
other corrections to scaling are small. 

We performed long spreading runs with a maximum 
time of 5 x 10 7 , system size 500 3 and vacancy concen- 
tration p = 0.5 for several infection rates A close to the 
phase transition. The resulting plot of N s vs P s is shown 
in Figure [SJ The data are averages over 10 5 to 10 6 dis- 



FIG. 6. (Color online) Effective critical exponent —<d/S — 
(d\nN s )/(d\nP s ) vs. P^ 2 calculated from the data in Fig. 
[5] The critical curve is marked by dots, and the dashed line 
is a linear extrapolation to P s = 0. 



order configurations with 1000 trials starting from ran- 
dom seed sites for each configuration. The figure shows 
that the relation between N s and P s indeed approaches 
a power law in the long-time (small P s ) limit. The fig- 
ure also indicates that the crossover to the asymptotic 
behavior is very slow. The asymptotic power law is only 
reached when P s falls well below 10 -3 which corresponds 
to times larger than 10 4 , implying that long simulations 
are required to determine the critical behavior. More- 
over, the mean-square radius of the active cluster at the 
crossover time is approximately 25, implying a total di- 
ameter of about 100. This means that simulations of 
systems with less than 100 3 sites will never reach the 
asymptotic critical behavior. 



D. Critical exponents 

To find the critical infection rate A c and to measure the 
finite-size scaling exponent /3/v± (the scale dimension of 
the order parameter), we define the effective (running) 
exponent (din N s ) / (din P s ) = —Q/S. It is related to the 
finite-size scaling exponent via Q/S = 3v±//3 — 2 [see 
cqs. (HJU) and (|2H)) ]. To avoid the uncertainties stemming 
from the unknown microscopic time scale to, we extrapo- 
late this effective exponent to P s = rather than t = oo. 
Figure O shows the effective exponent as a function of Pf 
with y — 1/2. (The value 1/2 was again chosen empiri- 
cally to permit an approximately linear extrapolation of 
the data with P s < 10~ 3 . Because of the slow crossover 
to the asymptotic regime, the value of y is much more 
uncertain than that of y in the clean case, see Fig. [2]) 

From Fig. [6l we conclude that A c = 2.6906(3) for a va- 
cancy concentration of p = 0.5. Extrapolating the effec- 
tive exponent to P s — yields Q/S = —0.42(3) where the 
error estimate largely stems from the uncertainty in A c 
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FIG. 7. (Color online) Effective critical exponents 5 = 
-(dlnP s )/(dlnln(t/t )) and 9 = (dlnN s )/(dlnln(t/t )) vs. 

1 /2 

P s for p = 0.5 and tp = 1. The critical curves are marked by 
dots, and the dashed lines are linear extrapolations to P a — 0. 



Value This work Ref. [54] 





1.90(4) 


1.84(2) 




0.98(6) 


0.99(2) 


p 


1.87(7) 


1.82(4) 


Df 


1.10(4) 


1.16(2) 


i> 


0.38(3) 


0.46(2) 




0.37(4) 


0.45(3) 


5 


5.0(2) 


4.0(2) 


e 


-2.1(2) 


-1.5(1) 



TABLE II. Critical exponents of the disordered three- 
dimensional contact process compared to results of the SDRG 
calculation [5J]. The upright numbers are directly given in 
Ref. [54|, the italic ones were calculated using scaling rela- 
tions. The fractal dimension Df — 3 — P/v±. 



hampered by the rather large uncertainties in ip and to- 
A better estimate can be obtained by considering the 
dependence of (dhxP s )/(d\) and (dhxN B )/(dX) on P s 
which takes the form 



(and the related uncertainty in y.) The statistical error is 
much smaller. The resulting finite-size scaling exponent 
is f}/v± — 1.90(4). An estimate for this exponent can also 
be obtained from analyzing the dependence of N s on R 
in a similar fashion. The data show additional curvature 
(corrections to scaling), thus giving the less precise value 
P/v x = 1.85(15). 

We now apply the same type of analysis to the logarith- 
mic derivatives of P S) N s , and R with respect to \n(t/to) 
to determine the values of the exponents 5, 0, and ip. 
This requires a value for the microscopic time scale to. 
Since an incorrect to would produce additional correc- 
tions to scaling, we estimated its value by minimizing 
the time dependence (P s dependence) of the effective ex- 
ponents 8 and 0. This yields to « 1.0(4). Figure [7] shows 
the resulting effective exponents 8 and as a function of 

1/2 

P?' for vacancy concentration p = 0.5. Extrapolating 
the data at the critical infection rate A c = 2.6906 to 
P s = (i.e., t = oo) gives the values 6 — 5.0(2) and 
= —2.1(2). Again, the error estimate is dominated 
by the uncertainty in A c (and the resulting uncertainties 
in y and to). The tunneling exponent ip can be deter- 
mined by combining the value of /3/u± with either 5 or 
0. We find ip = 0.38(3). Alternatively, ip can be ob- 
tained from the dependence of R on \n(t/to). As these 
data show additional corrections to scaling, the extrapo- 
lation to t = oo is difficult and leads to the less precise 
estimate ip = 0.41(5). 

To find the critical exponents v± and j3, we now study 
the dependence of (dlnP s )/(dX) and (dhiN s )/(dX) on 
ln(t/to) according to eq. (|2"TT) . This yields the exponent 
combination l/v±ip = 2.7(3). Combining this with the 
value for ip, we obtain u± = 1.0(2). This analysis is 



d In P ' d In N . iy/J f . 

Extrapolating the effective exponents to P s = as before, 
we obtain the values 1//3 =0.53(2) and 0.55(3) from the 
P s and N s data, respectively. Our final estimate of the 
order parameter exponent is thus j3 = 1.87(7). Combined 
with the finite-size scaling exponent, this yields v\_ = 
0.98(6). 

In Table [Til we compare our estimates for the criti- 
cal exponents with results of a numerical SDRG calcula- 
tion [54[ of the random transverse-field Ising model with 
up to 128 3 sites. This model is expected to be in the 
same universality class as the disordered contact process. 
All static exponents (above the dividing line in the ta- 
ble) agree within their error bars (though just barely in 
the case of /3/v±). In contrast, the tunneling exponent 
ip and the other exponents characterizing the time de- 
pendencies (below the dividing line) do not agree. This 
suggests that the uncertainties in determining the micro- 
scopic time scale to and, correspondingly, the microscopic 
energy scale Slo of the SDRG calculation may be respon- 
sible for the disagreement because to and flo do not in- 
fluence the static exponents. We note, however, that our 
raw data do not seem to be compatible with the values 
for 5 and calculated from the results of Ref. [54j even if 
we allow to to vary (unless one assumes a crossover from 
our 6 — 5 to 6 — 4 and from = 2 to = 1.5 at times 
t > 10 s beyond the range of our simulations). This can 
be seen in Fig. [8] where we compare our data to the func- 
tions P s - ln(i/i )" 4 and N s - ln(t/t ) -1 ' 5 with ln(i ) 
= 0, 1, 2, 3, and 4. We will return to this question in the 
concluding section. 



In(t) 

FIG. 8. (Color online) Time dependence of P 3 and N s for 
impurity concentration p = 0.5 and several A close to the 
transition compared to the predictions of the numerical SDRG 
of Ref. [54|. The dashed lines represent the functions P 3 ~ 
ln(^o)" 4 and N s ~ rn(t/*o) _1 ' S with ln(t ) = 0, 1, 2, 3, 4 
(bottom to top) and arbitrary prefactor. 



E. Universality of the critical behavior 



1 °10 1 10 2 10 3 10 4 10 5 10 6 10 7 

t 

FIG. 9. (Color online) Upper panel: N s vs. P s at criticality 
for several vacancy concentrations p. Lower panel: P s vs. t at 
criticality, demonstrating the crossover from the clean to the 
dirty critical behavior. The dashed line represents a power 
law with the clean critical exponent 5 — 0.7366 and arbitrary 
prefactor. 



So far, all results on the disordered contact process 
were for a vacancy concentration p = 0.5. We now ad- 
dress the question of whether or not the critical behavior 
is universal, i.e., independent of the disorder strength. 
The SDRG underlying the infinite-randomness scenario 
becomes exact only for infinitely strong disorder (in- 
finitely broad disorder distributions). Therefore, it can- 
not decide the fate of a weakly disordered system. How- 
ever, Janssen's perturbative renormalization group [l8j . 
which is controlled for weak disorder, shows runaway flow 
towards large disorder strength. Furthermore, Hoyos [55] 
showed that within an improved SDRG scheme, the dis- 
order always increases under renormalization, even if it 
is weak initially. These arguments support a universal 
scenario in which the critical behavior is independent of 
the disorder strength. 

To study the question of universality numerically, we 
performed simulations for vacancy concentrations p = 
0.2, 0.3, 0.4, and 0.6 in addition to the value 0.5. Re- 
peating the complete analysis as discussed in the pre- 
vious subsections for all values of p would have been 
prohibitively expensive in terms of computer time. We 
therefore focused on finding the finite-size scaling expo- 
nent from N s vs. P s plots analogously to Fig. [5J using 
somewhat shorter runs. The maximum time was at least 
3 x 10 6 for all vacancy concentrations, and the data are 



averages over at least 10 7 trials using systems of 500 3 or 
720 3 sites. The resulting critical curves are presented in 
the upper panel of Fig. In the low-P, (long-time) limit, 
all curves appear to be parallel, implying that Q/S and 
with it the finite-size scaling exponent /3/v± takes the 
same value for all vacancy concentrations. The figure 
also suggests that the weak-disorder curves (in particu- 
lar p — 0.2) have not fully crossed over to the asymptotic 
critical behavior. This is confirmed in the lower panel 
of Fig. [9] which presents a log-log plot of P s vs. t at 
criticality. While the stronger-disorder curves (p = 0.5, 
0.6) start to deviate from the clean critical power law at 
t * 10 3 to 10 4 (in agreement with the estimate discussed 
at the end of Sec . ITVC)) , the p = . 2 curve deviates appre- 
ciably only after t « 10 6 . (Note that these long crossover 
times also imply huge system sizes to reach the asymp- 
totic regime. For p — 0.2, the mean-square radius of the 
active cloud at the crossover time of 10 6 is about 200.) 

Our simulations thus show no indications of nonuniver- 
sal, continuously varying critical exponents. However, we 
cannot rigorously exclude that the exponents change for 
very weak disorder, because the extremely large crossover 
times between the clean and the dirty critical behavior 
prevent us from reaching the asymptotic regime in these 
cases. 
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FIG. 10. (Color online) P s vs. t for several infection rates A 
in the Griffiths region Al? < A < A c (vacancy concentration 
p — 0.5). The dashed lines are power-law fits of the long- 
time behavior to (|25[) . Inset: Resulting Griffiths dynamical 
exponent z' as a function of A. The solid line is a fit to (|26[1 . 



FIG. 11. (Color online) P s vs ln(t) for p = p c = 0.6883920 and 
A = 6.0. The dashed line represents the expected logarithmic 
long-time decay, P s ~ [ln(t)]~ <s with 8 = 0.188 and arbitrary 
prefactor. 



F. Griffiths region 

In order to investigate the Griffiths region A° < A < A c , 
we have also performed detailed simulations for infection 
rates A below but close to the critical rate A c . Figure [TOl 
presents the resulting survival probability P s as a func- 
tion of time t for vacancy concentration p — 0.5. The 
data are averages over at least 10000 disorder configura- 
tions with 1000 trials per configuration. The system size 
is 300 3 sites. For all infection rates shown, the long-time 
decay of the survival probability obeys (over several or- 
ders of magnitude in P s and/or t) the nonuniversal power 
law predicted by the rare region arguments of Sec. MI CI 

The Griffiths dynamical exponent z' can be found by 
fitting the long-time decay to (1251) . The resulting val- 
ues, shown in the inset of Fig. [TU1 demonstrate that z' 
diverges as A approaches the critical value A c = 2.6906. 
Fitting z' to the expected power law (|2"6"|) gives a value 
for the combination v±ip. The fit is not of particularly 
high quality, but the resulting value, vj_tp = 0.42(6), is 
in reasonable agreement with the value determined at 
criticality. 



G. Contact process at the lattice percolation 
threshold 

This subsection is devoted to the nonequilibrium phase 
transition of the diluted contact process across the lattice 
percolation threshold p c . In the phase diagram shown in 
Fig. [TJ this transition is marked by the vertical line at p c 
between A^ 1 = and the multicritical point. 

The contact process on a diluted lattice close to the 
percolation threshold can be understood by combining 
classical percolation theory with the properties of the su- 



percritical contact process on finite-size clusters [37J, l38| ■ 
Although its behavior follows the activated scaling sce- 
nario described in subsection IIII B[ the critical expo- 
nents of the percolation transition differ from those of 
the generic transition discussed in the preceding sections. 
Interestingly, they are completely determined by the val- 
ues of the classical lattice percolation exponents /3 C and 
v c which are known numerically with high accuracy [56[ . 

According to Refs. [U HH, the order parameter ex- 
ponent ,3 and the correlation length exponent y± of 
the nonequilibrium phase transition are identical to the 
corresponding lattice exponents, /3 = f3 c = 0.417 and 
vj_ = v c = 0.875. The tunneling exponent ip is given 
by the fractal dimension D c of the critical lattice per- 
colation cluster, ip = D c = 3 — fi c jv c = 2.523. As a 
result, the critical exponent S takes a very small value, 
5 = P/(y ± <ip) = 0.188. 

We performed spreading simulation runs (starting from 
a single active site) at p c = 0.6883920 48] and A = 6.0 to 
a maximum time of 5 x 10 5 . Because of the small value 
of S, these simulations are particularly time consuming. 
The resulting survival probability P s (averaged over 230 
disorder configurations with 200 trials per configuration) 
is presented in Fig. [TTJ The data are in agreement with 
the qualitative predictions of Refs. [37l |38| : a rapid ini- 
tial decay towards a quasi-stationary state, followed by 
a slow logarithmic time dependence due to the succes- 
sive dying out of the contact process on larger and larger 
lattice percolation clusters. The figure shows that the 
long-time behavior of our data is compatible with the 
predicted exponent value 5 = 0.188. However, we found 
it impossible to calculate a precise value of 6 directly 
from the simulation data because the long-time decay is 
extremely slow. 
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V. SUMMARY AND CONCLUSIONS 

To summarize, we performed large-scale Monte Carlo 
simulations of the contact process on site-diluted cubic 
lattices. We determined the infection rate-dilution phase 
diagram. It features two different nonequilibrium phase 
transitions, (i) the generic transition that occurs for di- 
lutions below the percolation threshold of the lattice and 
is driven by the dynamic fluctuations of the contact pro- 
cess, and (ii) the transition across the percolation thresh- 
old which is driven by the lattice geometry. 

Our simulation results show that the generic transi- 
tion is controlled by an infinite-randomness critical point 
for all dilutions investigated. It gives rise to ultra- 
slow activated (exponential) dynamical scaling instead 
of the power-law dynamical scaling at conventional crit- 
ical points. The corresponding logarithmic time depen- 
dencies of various observables at criticality required long 
simulation times and thus a huge numerical effort (in to- 
tal about 100000 CPU days on the Pegasus cluster at 
Missouri S&T). We determined the complete critical be- 
havior of the generic transition and found it to be univer- 
sal, i.e., independent of the disorder strength (dilution). 

The critical exponents are listed in Table HIT to gether 
with results of a numerical SDRG calculation |54| of the 
three-dimensional random transverse-field Ising model 
which is predicted to be in the same universality class. 
We were able to calculate reasonably accurate estimates 
for the static critical exponents including the finite-size 
scaling exponent f3/v±, the order-parameter exponent /3 
and the spatial correlation length exponent v±_ . The cor- 
relation length exponent satisfies the inequality dv±_ > 2, 
as is expected in a disordered system [57[. Our values 
for the static exponents agree with the corresponding nu- 
merical SDRG results of Ref. [54| within the given errors. 
In contrast, our result for the tunneling exponent if) as 
well as 5 and O do not agree with the values quoted in 
Ref. [54| | . Estimates of -ipi 5 an d O depend sensitively 
on the microscopic time scale t$ or, correspondingly, on 
the microscopic energy scale ilo in the numerical SDRG 
calculation, while the static exponents are independent 
of it. This suggests that uncertainties in the value of to 
or £!o may be responsible for the disagreement. However, 
even if we allow the value of to to vary, our P s and N s 
data do not seem to be compatible with the values of 5 
and O predicted by Ref. [54 1. 

The differences between our results and the numeri- 
cal SDRG calculation could either imply a real differ- 
ence in universality class, or they could mean that one or 
both sets of results represent effective rather than true 
asymptotic exponents. A full resolution of this question 
will likely require much more extensive simulations to- 
gether with a careful analysis of finite-size effects and, 
in particular, of the microscopic time/energy scale in the 
infinite-randomness scenario. 

In addition to the generic transition, we briefly stud- 
ied the nonequilibrium transition across the lattice per- 
colation threshold. The Monte Carlo results support the 



predictions of the theory developed in Refs. [37|,l38|: The 
dynamical critical behavior is of activated type with crit- 
ical exponents that are combinations of the classical lat- 
tice percolation exponents. Our simulations also allowed 
us to find with reasonable accuracy the location of the 
multicritical point separating the generic transition from 
the percolation transition (see Fig. [T|). However, as the 
dynamics is expected to be even slower than that of the 
generic transition, finding the true multicritical behav- 
ior appears to be beyond our current computational re- 
sources. 

We also obtained high-precision estimates for the crit- 
ical behavior of the three-dimensional DP universality 
class by performing simulations of the clean contact pro- 
cess in three dimensions. We employed the reweighting 
technique proposed by Dickman [47| to save computer 
time. By using large lattices of up to 999 3 sites and long 
times of up to 5 x 10 4 , we were able to compute the crit- 
ical exponents with unprecedented accuracy (see Table 
[I]). As the dynamics of the clean contact process is much 
faster than that of the disordered contact process, this 
part of the work took only a small fraction (about 5000 
CPU days) of the overall computer time. 

Let us now put our results into broader perspective. 
Our results for the critical behavior of the disordered con- 
tact process in three dimensions are in agreement with 
a general classification [H| of phase transitions in 
quenched disordered systems according to the effective 
dimensionality d e s of the defects and the lower critical 
dimension d~ of the problem. (A) If c? c ff < d~ , the crit- 
ical point is of conventional power-law type and accom- 
panied by exponentially weak Griffiths singularities. In 
class (B), which contains systems with d e g = d~ , the 
critical behavior is controlled by an infinite-randomness 
fixed point with activated scaling, accompanied by strong 
power-law Griffiths singularities. (C) For d c g > d~ , 
the rare regions can undergo the phase transition inde- 
pendently from the bulk system. This leads to a de- 
struction of the sharp phase transition by smearing [59j j . 
For the contact process with vacancies (point defects), 
d c ff = d~ = leading to class B. In contrast, the contact 
process with extended (line or plane) defects belongs to 
class C [H \($. 

We conclude by noting that the exotic critical behavior 
of the disordered contact process (in one, two, and three 
dimensions) may be responsible for the striking absence 
of directed percolation scaling in at least some of the ex- 
periments [12| ■ In view of the increased experimental ac- 
tivities in the area of absorbing-state transitions [l3l - fl6| , 
we hope that our theoretical results will help guiding the 
data analysis in further experiments. However, it must 
be pointed out that the extremely slow dynamics and 
narrow critical region will prove to be a challenge for the 
verification of the activated scaling scenario not just in 
simulations but also in experiments. Finally, we empha- 
size that our results are of importance beyond absorbing 
state transitions. The strong-disorder renormalization 
group predicts our transition to belong to a broad univer- 
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sality class that also includes, e.g., the three-dimensional 
random transverse- field Ising model [28l . [29l I41JJ . Conse- 
quently, the critical behavior found here should be valid 
for other systems in this universality class as well. 
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